function ben_measureDistances(varargin)
% BEN_MEASUREDISTANCES  Measures cortical thickness and relative distance
% from pia of soma and all laminar borders (L1/L2, L3/L5a, L5a/L5b).
%
% Points should be mouse-selected in the following order:
% 1. soma (right-click)
% 2. tangent to pia (left-click, right-click)
% 3. L1/2 border (left-click)
% 4. L3/5a border (left-click)
% 5. L5a/5b border (left-click)
% 6. white matter (right-click)
%
% Data are saved to the workspace as structure 'ben_distanceData'
%
% This new version expects the user to first right-click on the soma, 
% then draw a line parallel to the pia, roughly tangent, 
% by left-clicking for 1st point, then right-clicking for 2nd
% point. A helper line through the soma, perpendicular to the pia, is then
% drawn and the user is expected to click the following 5 points along this
% line: pia, L1/2 border, L3/L5a border, L5a/L5b border, white matter.
%
% optional arguments:
%   1. imgFilePath   path of image file, for record-keeping
%   2. mode          'legacy' or 'compatibility'
%
% Ben Suter 2009-10-09 -- initial new version with helper line, etc.
% Ben Suter 2000-10-12 -- prints structure ready for copy-and-paste
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[sx, sy] = ben_getpts; % get soma
if numel(sx) ~= 1
    beep;
    disp('ben_measureDistances: Must first select just the soma!');
    return;
end

[px, py] = getline; % get pia parallel, roughly tangent
if numel(px) ~= 2
    beep;
    disp('ben_measureDistances: Must define line parallel to pia with exactly two points!');
    return;
end

p1 = [ px(1) py(1) ];
p2 = [ px(2) py(2) ];
s = [ sx sy ];

% Find point on pial parallel closest to soma
m = p1 + dot(p2-p1, s-p1)*(p2-p1)/dot(p2-p1, p2-p1);

% Draw a line through soma, orthogonal to pia, from bottom of screen to top
yl = ylim(); 
slope = (s(2)-m(2)) / (s(1)-m(1));
bottomx = s(1) + (yl(1)-s(2)) / slope;
bottomy = yl(1);
topx = s(1) - (s(2)-yl(2)) / slope;
topy = yl(2);
helperLine = line([bottomx, topx], [bottomy, topy], 'LineStyle', ':', 'Color', 'g');

[x, y] = ben_getpts; % get pia, borders, wm (5 points)
if numel(x) ~= 5
    beep
    disp('ben_measureDistances: Must select exactly 5 points. Type help ben_measureDistances for instructions.');
    return
end

% Get rid of the helper line, it is no longer needed
delete(helperLine);

% Position of landmarks
pia = [x(1) y(1)];
L1L2 = [x(2) y(2)];
L3L5a = [x(3) y(3)];
L5aL5b = [x(4) y(4)];
wm = [x(5) y(5)];

% An inline function to compute distance between two vectors
distance = @(a, b) sqrt(dot(b-a, b-a));

% Distance between landmarks (we measure distance relative to pia)
corticalThickness = distance(wm, pia);
fracSoma = distance(s, pia) / corticalThickness;
fracL1L2 = distance(L1L2, pia) / corticalThickness;
fracL3L5a = distance(L3L5a, pia) / corticalThickness;
fracL5aL5b = distance(L5aL5b, pia) / corticalThickness;

ben_distanceData.corticalThickness = corticalThickness;
ben_distanceData.fracSoma = fracSoma;
ben_distanceData.fracL1L2 = fracL1L2;
ben_distanceData.fracL3L5a = fracL3L5a;
ben_distanceData.fracL5aL5b = fracL5aL5b;
ben_distanceData.somaPosition = [sx sy];
ben_distanceData.angleFromVertical = - atan(1/slope) * 180/pi;

if size(varargin,2) > 0
    ben_distanceData.imageFile = varargin{1};
else
    ben_distanceData.imageFile = '';
end

if size(varargin,2) > 1 
    if strcmp(varargin{2}, 'legacy') 
        generateLegacyOutput(pia, L1L2, s, wm);
    elseif strcmp(varargin{2}, 'compatibility')        
        % For backwards compatibility
        ben_distanceData.yL2 = fracL1L2 * corticalThickness;
        ben_distanceData.ysoma = fracSoma * corticalThickness;
        ben_distanceData.yWM = corticalThickness;
        ben_distanceData.yfrac = fracSoma;
    end
end

assignin('base', 'ben_distanceData', ben_distanceData);

% %%%%%%% Optionally print results in .m-file ready format %%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
doDataPrint = true;
if doDataPrint
    prefix = 'map.image.distanceData.';
    disp('% ********** Image: distance measurement **********')
    fnames = fieldnames(ben_distanceData);
    for n = 1 : size(fnames,1)
        val = ben_distanceData.(fnames{n});
        if max(size(val)) == 1 % scalars
            disp([prefix fnames{n} ' = ' num2str(val) ';'])
        else % vectors
            if ischar(val)
                disp([prefix fnames{n} ' = ''' num2str(val) ''';'])
            else
                disp([prefix fnames{n} ' = [' num2str(val) '];'])
            end
        end
    end
end

end  % end of main function

function generateLegacyOutput(pia, L1L2, soma, wm)
% GENERATELEGACYOUTPUT Calculates distance measurements by the original
% method, from the positions of pia, L1/L2 border, soma, wm selected using
% the new method, and puts these into the workspace in the original
% (legacy) format.
%
% Calculates x, y, d, yL2, ysoma, yWM, yfrac and dumps these into the
% workspace.
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xy = [ pia; L1L2; soma; wm; ];
x = xy(:,1); y = xy(:,2);

% From here on this code is exactly the same as in the old measureDistances
% function

distanceData.x = x;
distanceData.y = y;

disp(' ');
disp('measureDistances: X positions of points:');
if max(abs(x))>1
    disp(num2str(round(x')));
else
    disp(num2str(x'));
end

disp(' ');
disp('measureDistances: Y positions of points:');
if max(abs(y))>1
    disp(num2str(round(y')));
else
    disp(num2str(y'));
end

if numel(x)<2
    beep
    disp(' ');
    disp('measureDistances: need to select at least 2 points to measure a distance.');
    disp(' ');
    return
end

for n = 2:size(x,1);
    dx = x(n) - x(n-1);
    dy = y(n) - y(n-1);
    d(n) = sqrt(dx^2 + dy^2);
end
distanceData.d = d;

disp(' ');
disp('measureDistances: distances between points:');
if max(abs(d))>1
    disp(num2str(round(d)));
else
    disp(num2str(d));
end

disp(' ');
disp('measureDistances: cumulative distances between points:');
if max(abs(d))>1
    disp(num2str(round(cumsum(d))));
else
    disp(num2str(cumsum(d)));
end
disp(' ');

% yfrac calculation
if numel(d) == 4
    if max(abs(d))>1
        D = cumsum(d);
        ysoma = D(end-1);
        yWM = D(end);
        yfrac = round(1000*(ysoma/yWM))/1000;
        disp('measureDistances: yfrac:');
        disp(num2str(yfrac));
        
        distanceData.yL2 = D(end-2);
        distanceData.ysoma = ysoma;
        distanceData.yWM = yWM;
        distanceData.yfrac = yfrac;
    else
        disp('measureDistances: cumulative distances between points:');
        disp(num2str(cumsum(d)));
    end
    disp(' ');
end

assignin('base', 'distanceData', distanceData);
end
